Geodynamically corrected Pliocene shoreline elevations in Australia consistent with midrange projections of Antarctic ice loss

The Mid-Pliocene represents the most recent interval in Earth history with climatic conditions similar to those expected in the coming decades. Mid-Pliocene sea level estimates therefore provide important constraints on projections of future ice sheet behavior and sea level change but differ by tens of meters due to local distortion of paleoshorelines caused by mantle dynamics. We combine an Australian sea level marker compilation with geodynamic simulations and probabilistic inversions to quantify and remove these post-Pliocene vertical motions at continental scale. Dynamic topography accounts for most of the observed sea level marker deflection, and correcting for this effect and glacial isostatic adjustment yields a Mid-Pliocene global mean sea level of +16.0 (+10.4 to +21.5) m (50th/16th to 84th percentiles). Recalibration of recent high-end sea level projections using this revised estimate implies a more stable Antarctic Ice Sheet under future warming scenarios, consistent with midrange forecasts of sea level rise that do not incorporate a marine ice cliff instability.


Supplementary Text
Figs. S1 to S11 Table S1 References

Supplementary Text S1 Optimal models of present-day dynamic topography
Prior to investigating the impact of time-evolving dynamic topography on Australian palaeo sea-level data, we must first construct Earth models that satisfy all available geodynamic observables that relate to the present-day condition of the mantle.To do so, we adopted the approach of Richards et al. (34), in which mantle density structures are sought that optimise global fits to independent observations of dynamic topography, geoid height anomalies, and core-mantle boundary (CMB) excess ellipticity.
These density models are constructed in a near-identical manner to the input temperature fields used in our time-dependent simulations (see Section M2 above).Density above 400 km is derived from a modified version of the RHGW20 temperature and density model (61).Deeper than 300 km and outside dense basal sections of LLVPs, the same anelastically corrected and smoothed pyrolite look-up table is used to convert seismic velocities from five different V S tomographic models [LLNL-G3D-JPS (37); S40RTS (38); SAVANI (39); SEMUCB-WM1 (40); TX2011 (41)] into temperature and density.
Mid-mantle V S structure is also high-pass filtered over the 1000-2000 km depth range and the two different density parameterisations are smoothly merged between 300 km and 400 km depth by taking their weighted average.
The principal difference between the density inputs used in our instantaneous and time-dependent mantle flow models is that, rather than setting a fixed density jump for compositionally distinct material within the basal LLVP layer, a range of non-pyrolitic compositional endmembers are explored.These include mid-ocean ridge basalt, chondrite-enriched basalt, and iron-enriched pyrolite [see Table 1 in Richards et al. (34)].In cases where the composition is intermediate between pyrolite Table S1: Optimal LLVP basal layer parameters for each combination of mantle density model and radial viscosity profile.Layer composition is assumed to be a mechanical mixture of pyrolite (89) and the quoted chemical component: MORB = mid-ocean ridge basalt (89); CEB = chondrite-enriched basalt (90); FSP = iron-enriched pyrolite (91,92).

Density
Viscosity Thickness (km)  and a given endmember, properties appropriate for a mechanical mixture of the two components are calculated using the Voigt-Reuss-Hill approximation to average the elastic moduli and generate an appropriate anelastically corrected and smoothed V S lookup table.For all models where LLVP basal layer composition is considered distinct from ambient mantle (i.e., the mantle is not 100% pyrolite), temperatures and densities are determined separately for the two components and then combined into a single array.For each candidate compositional endmember, density models are generated for enrichments of [0, 10, ..., 100]% and for dense basal layer thicknesses of 50 km and 100-900 km in 100 km increments.In each case, the mean density anomaly of the combined intra-and extra-LLVP regions is normalised to zero such that the mean density remains equal to PREM (93).By constructing these density structures from each of our five deep mantle V S models, we generate an ensemble of 505 mantle density inputs.
For each density input, we predict geoid undulations, surface dynamic topography, and CMB dynamic topography up to a maximum spherical harmonic degree, l max = 30, using the instantaneous flow kernel methodology outlined in Corrieu et al. (94).This formulation accounts for the effects of both compressibility and self-gravitation [see Appendix B of Richards et al. (34) for details].In each case, we assess the sensitivity of our results to Earth's uncertain viscosity structure by computing instantaneous mantle flow for three previously published radial profiles that are constrained by geoid, heat flow and glacial isostatic adjustment observations: S10 (35); F10V1; and F10V2 (36).For each pair of tomographic and radial viscosity inputs (15 in total), we then determine the optimal LLVP basal layer thickness and composition that fits observed geoid undulations, surface dynamic topography, and CMB excess ellipticity.Successful models must also exhibit an intrinsic density anomaly ≥ 2% for chemically distinct material within the LLVP, which is the estimated lower bound for long-term preservation of chemical heterogeneity (95,96,97) [see Richards et al. (34) for misfit parameterisation].
The 15 best-fitting mantle density models obtained via this optimisation procedure yield LLVP basal layer thicknesses between 50 km and 210 km and 2.0-3.9%intrinsic compositional density differences.There is evidence for a negative trade-off between these two LLVP parameters.Nevertheless, for the purposes of this study, it is important to note that their predicted present-day dynamic topography fields remain consistent (1σ ∼190 m; Figure S1) and that they yield especially good agreement with oceanic residual depth anomalies around the margins of Australia (regional misfit is 25-65% smaller than the global average).These characteristics suggest that our present-day mantle density models are both robust and mutually compatible beneath this region, thereby validating their use in the time-dependent simulations that are central to this study.

S2 Modelling time-dependent mantle flow
Having verified that our optimised mantle density models yield good agreement between predicted and observed present-day dynamic topography around Australia, we next incorporate these inputs into time-dependent simulations to hindcast (or 'retrodict') changes in dynamic topography since Mid-Pliocene times.We chose to implement these retrodictions using the convection code ASPECT, which required us to modify the V S -to-density conversion methodology used in our instantaneous models to allow for reasonable computation times.First, rather than inferring lower mantle density from composition and temperature using Perple X lookup tables, we instead apply a radial thermal expansivity profile.Our profile is determined by fitting a three-part linear model to laterally averaged depth-dependent expansivities obtained from our hybrid V S -to-temperature conversion [i.e., the Richards et al. (61) anelastic parameterisation shallower than 400 km, the Perple X-derived pyrolite lookup table deeper than 300 km, and a weighted combination of the two between 300 km and 400 km depth; Figure S2].Secondly, our ASPECT simulations are carried out under the Boussinesq approximation (i.e., they are incompressible).We therefore adjust our models by removing depth-dependent increases in temperature and density caused by adiabatic compression and by reducing the intrinsic compositional density jump for anomalous LLVP material to 0-132 kg m −3 (i.e., 0-4% of the incompressible reference density, ρ 0 = 3300 kg m −3 ).This adjustment minimises the difference between compressible and incompressible model predictions of dynamic topography, with compressible models yielding 10-15% higher dynamic topography at long wavelengths (l = 1-5) and minimal offset at shorter wavelengths.These differences are negligible in the context of this study, since changes in dynamic topography since 3 Ma are dominated by short-wavelength signals driven by shallow mantle dynamics.Although compressibility has a larger effect on the predicted non-hydrostatic geoid, predicted rates of change exhibit minimal differences and, since radial surface displacement has a stronger impact on relative sea-level variation than perturbations to the equipotential, the impact on our sea-level predictions is modest (generally < 0.1 m).
Despite details of the density fields exhibiting substantial differences (in keeping with the variations between tomographic models that underpin them), we find that, in all cases, the long-wavelength pattern of mantle flow is dominated by near-vertical lower mantle downwelling beneath Australia (Figures S3-S4).Similarly, shallow mantle upwellings are predicted to occur beneath Cape Range, Tasmania, the South Eastern Highlands, and Cape York in all models; however, the connection between these features and deep mantle structure is not clear in all cases (e.g., compare Figures S3d-e and Figures S4d-e   Some systematic differences do, however, emerge in the long-wavelength pattern of dynamic to- Changes in the thickness and intrinsic density of the basal LLVP layer have a minor but detectable impact on predicted dynamic topography patterns, with thicker, denser endmembers generally producing less uplift in the southwest and more in the northeast (Figures S6 and S7).This pattern is to be expected since the southwestern edge of the Pacific LLVP strikes approximately northwest-southeast beneath Cape York, with variations in intrinsic LLVP density therefore influencing the gradient in dynamic topography perpendicular to this boundary.Modifying the plate motion model from MORVEL   to S12 has a larger impact on predicted dynamic topography change, with the latter producing higher rates of uplift and subsidence (Figures S6 and S7).Altering the plate motion model also induces substantial changes in the pattern of dynamic topography evolution; however, these differences depend strongly on the assumed viscosity and density model, reflecting associated changes in the direction and rate of flow in the shallow mantle.
Density models derived from S40RTS, LLNL-G3D-JPS, and SEMUCB-WM1 yield the best agreement with observed MPWP sea-level marker elevations, while the optimal radial viscosity profile covaries with input tomography and assumed LLVP structure.Best-fitting basal LLVP layer parameters combine relatively modest intrinsic density anomalies of 1-2% with upper bound layer thicknesses of ∼200 km.In all cases, significantly improved fit between predicted and observed sea-level marker elevations is obtained for the MORVEL plate motion model in comparison to S12 (r is on average ∼0.2 higher; Figures S6 and S7).

S3 Modelling Pliocene-to-Recent GIA-induced sea-level change
Australia is in the far field of both the former Fennoscandian and Laurentian ice sheets, the presentday Greenland and West Antarctic Ice Sheets, and most marine-based sectors of the East Antarctic Ice Sheet.Thus, ongoing continental levering in response to the last deglaciation is expected to be the dominant process responsible for introducing geographic variability into GIA-related relative sea-level change around the continent.There are two main sources of uncertainty in our GIA models that must be considered -the ice history and the mantle viscosity structure.
The configuration of ice sheets at the end of the MPWP is uncertain and different assumptions introduce changes to the amplitude of predicted relative sea-level variations (Figure S8).In particular, ice history 'B' (MPWP ice configuration equivalent to that of the present day minus the Greenland and West Antarctic ice sheets) yields predicted sea-level marker elevations that are up to ∼70% more negative than those obtained with ice history 'A' (MPWP ice configuration identical to the present day).This difference arises due to additional continental levering associated with the rapid increase in ice volume at 2.95 Ma in history A, but has minimal impact on our compilation of sea-level markers since they are mostly located close to the coastlines where the net contribution of this effect is close to zero.Differences in the assumed rheological structure of the mantle cause subtle modifications to the predicted pattern of marker elevation change (Figure S8).Switching from the 120p55 radial profile to 80CV2V2 increases the amplitude of both positive and negative deviations by up to 50%.This result reflects both an increase in levering amplitude associated with the thinner lithosphere in 80CV2V2 and an overall stronger GIA response to deglacial ice-mass change caused by the reduced lower mantle viscosity in this model.
Nonetheless, these models indicate that no plausible combination of ice sheet configuration or mantle rheological structure can be responsible for the observed spatial variability in MPWP sealevel marker elevations around Australia.The predicted amplitudes of elevation change are almost two orders of magnitude too small and the average correlation with the marker elevations is only r ∼ −0.04.Furthermore, the effect of introducing lateral viscosity variations into GIA predictions has been shown to be relatively minor around Australia (67).This factor suggests that our use of simple, purely radial Earth structure to explore GIA uncertainties is sufficient.Nevertheless, despite its limited ability to account for spatial variability in sea-level indicator elevations, correcting for GIA processes does increase our final MPWP GMSL estimate by ∼4 m.
Although we do not account for East Antarctic ice loss in the full GIA calculations outlined above, we have quantified the potential impact of melting certain marine-based sectors on measured RSL at our MPWP sea-level marker sites.These tests simulate steady loss of ice from Prydz Bay (∼1.5 m GMSLE) and Wilkes Subglacial Basin (∼3 m GMSLE) over the course of 5 kyr using the same modelling framework as the full simulations and the 120p55 viscosity profile.We find that RSL deviates from GMSL by less than a metre at all five sites.This result indicates that our conclusions are not materially impacted by our omission of spatially variable sea-level signals that could be caused by Mid-Pliocene ice loss from different sectors of the East Antarctic Ice Sheet.

S4 Bayesian Gaussian process-based inversion for MPWP GMSL
By framing estimation of MPWP GMSL as a Bayesian Gaussian process-based regression problem, we are able to constrain the posterior probability distribution of its value given: i) prior knowledge of its likely magnitude and temporal variation; ii) constraints on the age and palaeo-water depth of individual sea-level markers; and iii) our ensemble of GIA and dynamic topography predictions.This inference process permits GMSL solutions that produce temporal evolution consistent with the timescale and variance parameters of the underlying Gaussian process, while, at any individual timestep, model realisations that minimise spatial variability in sea level are naturally favoured.Consequently, the resulting posterior distributions of individual input parameters yield useful information concerning which density and viscosity models best account for MPWP sea-level indicator elevations.They therefore provide new constraints on Earth's internal structure as well as some quantification of related uncertainties in estimated GMSL.
The posterior distribution for GIA-induced elevation changes closely resembles the prior in terms of both median value and uncertainty range.This result is to be expected since the elevation changes are much smaller than those associated with both the dynamic topography correction and the observational data spread (Figure S9).GIA-induced elevation changes are least well-constrained in the interior of the continent (Figure S9d).The maximum a posteriori (MAP) estimates are ∼ 0.1 for the GIA viscosity profile index and ∼ 11 m GMSLE for the ice history, indicating a slight preference for the 120p55 viscosity model and an MPWP ice sheet configuration in which large portions of the Greenland and West Antarctic ice sheets were initially absent (Figures S10l and S10m).Nevertheless, uncertainties on these two parameters are sufficiently large that no firm conclusions can be drawn from this analysis.In contrast, the posterior distribution for predicted elevation change due to dynamic topography is significantly different to the prior distribution.The median dynamic topography change is closest to that expected for a model based on the S40RTS tomography model, the F10V1 radial viscosity profile, and the MORVEL plate motion model, with an LLVP dense basal layer thickness of ∼ 100 km and intrinsic density contrast of ∆ρ C ∼ +2.0% (Figure S9).Spatial variations in elevation change uncertainty are variable, but are particularly pronounced in southeastern Australia and close to the Australian-Antarctic Discordance where there are substantial differences in the shallow to mid-mantle V S structure of the input tomographic models (Figure S9e).Dynamic topography dominates over GIA when the two fields are combined, as demonstrated by the similarity of panels b versus c, and e versus f in Figure S9.Posterior distributions for individual dynamic topography parameters and parameter pairs suggest that observed sea-level marker elevations are most consistent with mantle flow models based predominantly on S40RTS (MAP index value ∼ 0.3), with modest contributions from LLNL-G3D-JPS and/or SAVANI (MAP values ∼ 0.2 in both cases; Figures S10a-e and S11a-j and ca).Simulations with radial viscosity parameterisations similar to F10V1 are clearly optimal (MAP index value ∼ 0.5); however, composites of F10V1 and the other profiles investigated here, F10V2 and S10, also yield high posterior  S10f-h and S11u, aa, ab, and cb).The index for the plate motion model has the most tightly constrained posterior distribution of all input parameters, further demonstrating that the MORVEL plate motion model produces better agreement with our observations than the S12 model (Figures S10i and S11ac-aj, as, bb, bi, and bw).The posterior distributions for the two parameters controlling the properties of the dense basal layer within the LLVP (i.e., layer thickness and ∆ρ C ) are less well constrained and exhibit some negative trade-off, with thick, less dense layers yielding similar posterior probability to thin, denser layers (Figure S11bc).These distributions are also skewed towards the lower bound of each parameter, with MAP estimates of ∼ 25 km for layer thickness and ∼ +1.5% for ∆ρ C , but median estimates of ∼ 100 km and ∼ +2%, respectively (Figure S10j-k).

probability (Figures
It is important to note that parameters controlling the dynamic topography predictions do not trade off with those underpinning the GIA calculations (Figure S11).This result suggests that, in future work, it may be possible to disentangle the two contributions even when they have more comparable magnitudes, for example, in near-field locations or over shorter timescales (e.g., since the Last Interglacial).Although most of the dynamic topography parameters have well constrained posteriors, larger uncertainties for LLVP-related values likely reflect limited impact of the Pacific LLVP on mantle flow patterns beneath Australia.These parameters might be better determined by applying our analysis in other regions, such as Southern Africa, where LLVP properties will have a much stronger effect on the rate and pattern of dynamic topography change.
χG = global geodynamic misfit [Equation 7 of Richards et al. (34)]; r Aus A = Pearson's correlation coefficient for present-day dynamic topography and residual depth anomalies around Australian margins (see transect in Figure 1 of the main text and Figure S1); χ Aus A = misfit between present-day dynamic topography and residual depth anomalies around Australian margins [Equation 5 of Richards et al. (34)].Models in bold = those shown in Figures 1, S3 and S4; models in italics = those shown in Figure S1.

Figure S1 :
Figure S1: Consistency between predicted present-day dynamic topography and measured oceanic residual depth anomalies.(a) Predicted dynamic topography (water-loaded offshore, air-loaded onshore) for mantle density structure derived from SEMUCB-WM1 tomographic model (40) and the F10V2 mantle viscosity profile (36) expanded up to spherical harmonic degree, lmax = 30.Coloured circles/triangles = spot measurements of oceanic residual depth (86) (a common proxy for observed dynamic topography at present day); thick black line = location of transect shown in panel (e).(b) Same for LLNL-G3D-JPS tomographic model (37) and S10 mantle viscosity profile (35).(c) Mean dynamic topography across all 15 combinations of tomography and viscosity model.(d) Dynamic topography uncertainty (1σ) across all 15 models.(e) Predicted versus observed present-day dynamic topography along NW-to-SE transect.Red line/band = mean prediction with 1σ uncertainties; blue line = individual result for SEMUCB-WM1 and F10V2 combination; purple line = same for LLNL-G3D-JPS and S10; circles/triangles with error bars = spot measurements of residual depth and uncertainties (86); grey line = spherical harmonic fit to spot measurements (lmax = 30) with grey band representing range of values within a 500 km-wide swath perpendicular to transect.
beneath Cape Range and Tasmania).Despite substantial discrepancies in shortwavelength mid-mantle buoyancy structure, the pattern of predicted Mid-Pliocene-to-Recent dynamic topography change is consistent across our model ensemble, highlighting the dominant influence of the shallow mantle and the relative insensitivity of the vertical surface motions to deeper structure except at the longest wavelengths (Figures S5-S7).

Figure S3 :
Figure S3: Cross-sections through convection simulation based on SEMUCB-WM1 tomographic model (40) and F10V1 viscosity profile (36).(a) Grid = radial component of present-day mantle flow at 300 km depth; arrows = tangential component of flow; blue lines = plate boundaries.(b) Same as panel (a), except at 2700 km.(c-h) Cross-sections along transects shown in panels (a-b).Temperature anomaly given by red-blue colourscale; arrows = flow velocity in plane of transect; grey region = lithosphere (defined using depth to 1200 • C isotherm).

Figure S4 :
Figure S4: Cross-sections through convection simulation based on LLNL-G3D-JPS tomographic model (37) and F10V2 viscosity profile (36).Panels follow Figure S3.pography change predicted by the different models.SEMUCB-WM1-based models feature stronger upwelling along the western Australian margin, leading to greater uplift at Cape Range and in the Perth Basin, while LLNL-G3D-JPSand S40RTS-based models produce more pronounced uplift over the South Eastern Highlands (Figure S5).Models adopting the F10V1 viscosity profile yield the most uplift along the Eastern Highlands and relatively modest uplift along the western margin, with the reverse being true for S10-based counterparts.

Figure S5 :
Figure S5: Predicted Pliocene-to-Recent dynamic topography change as a function of tomography and viscosity input combination.(a) Predicted change in elevation of Mid-Pliocene Warm Period sea-level markers for LLNL-G3D-JPS tomographic model and S10 viscosity profile.Circles = Mid-Pliocene median uncorrected GMSL estimates (i.e., present-day elevation + palaeo-water depth; Table 2 in main text); r = correlation coefficient between predicted dynamic topography change and uncorrected GMSL estimates.(b) Same for LLNL-G3D-JPS and F10V1.(c) Same for LLNL-G3D-JPS and F10V2.(d-f) Same for TX2011.(g-i) Same for SEMUCB-WM1.(j-l) Same for S40RTS.(m-o) Same for SAVANI.Note that all models include a 100 km-thick, basal LLVP layer with ∆ρc = +2% intrinsic density contrast.Plate motion is set by the MORVEL plate velocity model (77).

Figure S7 :
Figure S7: Predicted Pliocene-to-Recent dynamic topography change as a function of intrinsic density contrast of basal LLVP layer and plate velocity model for the LLNL-G3D-JPS tomographic model and F10V2 viscosity pairing.Panels follow Figure S6.

Figure S8 :
Figure S8: Predicted GIA-induced change in MPWP sea-level marker elevations as a function of assumed mantle viscosity and ice history.(a) Predicted change in elevation of MPWP sea-level markers for 120p55 viscosity profile and ice history A, in which MPWP ice volume is equivalent to that at the present day.(b) Same for 120p55 viscosity profile and ice history B, in which Greenland and West Antarctic ice sheets are absent during the MPWP (∼+14 m GMSLE compared to present-day).(c-d) Same for 80CV2V2 viscosity profile.

Figure S9 :
Figure S9: Comparison of posterior GIA and dynamic topography elevation change predictions.(a) Posterior median for predicted GIA.Circles = Mid-Pliocene median uncorrected GMSL estimates (i.e., present-day elevation + palaeo-water depth; Table 2 in main text).(b) Same for dynamic topography.(c) Same for combined prediction.(d) Posterior uncertainty for predicted GIA (i.e., half of the difference between 16 th and 84 th percentile values).Circles = uncertainty in Mid-Pliocene GMSL estimates (i.e., combined elevation and palaeo-water depth uncertainty; Table 2 in main text).(e) Same for dynamic topography.(f) Same for combined prediction.Note change in colour bar ranges between GIA and dynamic topography panels due to substantially lower amplitudes in the former.

Figure S11 :
Figure S11: Trade-offs between model parameters.(a-bz) Joint posterior probability distributions (kernel density estimates) of model parameters.Solid contour = 84 th percentile; dashed contour = median value; dotted black lines = 16 th percentile; red circle = maximum a posteriori (MAP) value.Gray shading = region of parameter space that cannot be sampled since neither sum of tomographic model indices nor sum of dynamic topography viscosity profile indices can exceed unity.(ca) Posterior sample density (normalised so maximum value is unity) as a function of three tomographic model indices with highest MAP values (S40RTS, SAVANI, and LLNL-G3D-JPS).(cb) Same for indices of viscosity profiles used in dynamic topographic predictions.